3D Relativistic Magnetohydrodynamic Simulations of Current-Driven 
Instability. I. Instability of a static column 



Yosuke Mizuno 1,2 , Yuri Lyubarsky 3 , Ken-Ichi Nishikawa 1,2 , and Philip E. Hardee 4 

ABSTRACT 

We have investigated the development of current-driven (CD) kink instability 
through three-dimensional relativistic MHD simulations. A static force- free equilib- 
rium helical magnetic configuration is considered in order to study the influence of the 
initial configuration on the linear and nonlinear evolution of the instability. We found 
that the initial configuration is strongly distorted but not disrupted by the kink insta- 
bility. The instability develops as predicted by linear theory. In the non-linear regime 
the kink amplitude continues to increase up to the terminal simulation time, albeit at 
different rates, for all but one simulation. The growth rate and nonlinear evolution of 
the CD kink instability depends moderately on the density profile and strongly on the 
magnetic pitch profile. The growth rate of the kink mode is reduced in the linear regime 
by an increase in the magnetic pitch with radius and the non-linear regime is reached 
at a later time than for constant helical pitch. On the other hand, the growth rate 
of the kink mode is increased in the linear regime by a decrease in the magnetic pitch 
with radius and reaches the non-linear regime sooner than the case with constant mag- 
netic pitch. Kink amplitude growth in the non-linear regime for decreasing magnetic 
pitch leads to a slender helically twisted column wrapped by magnetic field. On the 
other hand, kink amplitude growth in the non- linear regime nearly ceases for increasing 
magnetic pitch. 

Subject headings: instabilities - MHD - methods: numerical - galaxies: jets 

1. Introduction 

Relativistic jets occur in active galactic nuclei (AGN) (e.g., Urry & Padovani 1995; Ferrari 
1998; Meier et al. 2001), occur in microquasars (e.g., Mirabel & Rodriguez 1999), and are thought 
responsible for the gamma-ray bursts (GRB) (e.g., Zhang & Meszaros 2004; Piran 2005; Meszaros 
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2006). The most promising mechanism for producing relativistic jets involves magnetohydrody- 
namic acceleration from an accretion disk around a black hole (e.g., Blandford 1976; Lovelace 1976; 
Blandford & Payne 1982; Fukue 1990; Meier 2005; Narayan et al. 2007), and/or involves the 
extraction of energy from a rotating black hole (Penrose 1969; Blandford & Znajek 1977). 

General relativistic magnetohydrodynamic (GRMHD) codes have been used to study the 
extraction of rotational energy from a spinning black hole, i.e., Blandford-Znajek mechanism (Koide 
2003; Komissarov 2005; McKinney 2005; Komissarov & McKinney 2007; Komissarov & Barkov 
2009) and from an accretion disk, i.e., Blandford-Payne mechanism (Blandford & Payne 1982). 
GRMHD codes also probe the formation of GRB jets in collapsars (e.g., Mizuno et al. 2004a, 2004b; 
Liu et al. 2007; Barkov & Komissarov 2008; Stephens et al. 2008; Nagataki 2009; Komissarov & 
Barkov 2009) and are used to study the propagation of jets injected from an unresolved Keplerian 
disk (e.g., Komissarov et al. 2007, 2008; Tchekhovskoy et al. 2008, 2009). GRMHD simulations 
with a spinning black hole indicate jet production with a magnetically dominated high Lorentz 
factor spine with v ~ c, and a matter dominated sheath with v > c/2 possibly embedded in 
a lower speed, v <C c, disk/coronal wind (e.g., De Villiers et al. 2003, 2005; Hawley & Krolik 
2006; McKinney & Gammie 2004; McKinney 2006; McKinney & Narayan 2007a, 2007b; McKinney 
& Brandford 2008; Beckwith et al. 2008; Hardee et al. 2007). Circumstantial evidence such as 
the requirement for large Lorentz factors suggested by the TeV BL Lacs when contrasted with 
much slower observed motions (Ghisellini et al. 2005) suggests such a spine-sheath morphology, 
although alternative interpretations are also possible (Georganopoulos & Kazanas 2003; Bromberg 
& Levinson 2008; Stern & Poutanen 2008; Giannios et al. 2009). 

Numerical simulations indicate that the jet spine and sheath are accelerated and, in part, 
collimated by strong magnetic fields twisted in the rotating black hole ergosphere and in the ac- 
cretion disk, respectively. The helically twisted magnetic fields are expected to be dominated by 
the toroidal component in the far zone as a result of jet expansion because the poloidal component 
falls off faster with expansion and distance. In configurations with strong toroidal magnetic field 
the current driven (CD) kink mode is unstable. This instability excites large-scale helical motions 
that can strongly distort or even disrupt the system. For static cylindrical force-free equilibria, 
the well-known Kruskal-Shafranov criterion states that the instability develops if the length of the 
column, £, is long enough for the field lines to go around the cylinder at least once (e.g. Bateman 
1978): \B P /B$\ < 1/2ttR. This criterion suggests that jets are unstable beyond the Alfven surface 
because of B p < B^ and £ > Rj at the Alfven surface, where Rj is the cylindrical radius of the 
jet. However, rotation and shear motions could significantly affect the instability criterion. For 
relativistic force- free configurations, the linear instability criteria were studied by Istomin & Pariev 
(1994, 1996), Begelman (1998), Lyubarskii (1999), Tomimatsu et al. (2001) and Narayan et al. 
(2009). 

The linear mode analysis provides conditions for the instability but says little about the impact 
the instability has on the system. The instability of the potentially disruptive kink mode found 
from a linear analysis must be followed into the non-linear regime. The nonlinear development 
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of the CD kink instability is different for external and internal kink modes (e.g., Bateman 1978). 
External disruptive instability develops in a current carrying plasma column surrounded by a 
vacuum magnetic field. In astrophysical systems, some amount of plasma is present everywhere. In 
this case an internal kink instability develops inside a resonant surface on which the helicity of the 
magnetic field lines matches that of the eigenmode. In dissipationless MHD, the internal instability 
is known to saturate. Finite resistivity makes the internal kink mode disruptive but the time-scale 
significantly exceeds the Alfven time-scale and may be too long to disrupt a fast outflow. In this 
case helical structures may be formed in the flow. Recently helical structures have been found in 
non-relativistic/relativistic simulations of magnetized jets (e.g., Lery et al. 2000; Ouyed et al. 2003; 
Nakamura & Meier 2004; Nakamura et al. 2007; Moll et al. 2008; McKinney & Blandford 2008; 
Carey & Sovinec 2009). 

Twisted structures are observed in many AGN jets on sub-parsec, parsec and kiloparsec scales 
(e.g., Gomez et al. 2001; Lobanov & Zensus 2001). Dissipation processes may result in relaxation 
to helical equilibrium as evoked by Konigl & Choudhuri (1985) for helical structures in AGN jets. 
In the absence of CD kink instability and resistive relaxation, helical structures may be attributed 
to the Kelvin-Helmholtz (KH) helical instability driven by velocity shear at the boundary between 
the jet and the surrounding medium (e.g., Hardee 2004, 2007) or triggered by precession of the jet 
ejection axis (Begelman et al. 1980). It is still not clear whether current driven, velocity shear driven 
or jet precession is responsible for the observed structures, or whether these different processes are 
responsible for the observed twisted structures at different spatial scales. 

This is the first in a series of papers in which we study kink instability in relativistic systems. 
By relativistic we mean not only relativistically moving systems but any with magnetic energy 
density comparable to or greater than the plasma energy density, including the rest mass energy. 
In this paper, we present 3D results of the CD kink instability of a static plasma column. We 
start from static configurations because in the case of interest, the source of the free energy is the 
magnetic field, not the kinetic energy. Therefore static configurations (or more generally rigidly 
moving flows considered in the proper reference frame) are the simplest ones for studying the basic 
properties of the kink instability. At the next stage, we will investigate the influence of shear 
motions and rotation on the stability and nonlinear behavior of jets. This article is organized as 
follows. We describe the numerical method and setup used for our simulations in §2, present our 
results in §3, discuss the astrophysical implications in §4, and present some numerical tests in the 
Appendix. 



2. Numerical Method 

In order to study time evolution of the CD kink instability in the relativistic MHD (RMHD) 
regime, we use the 3D GRMHD code "RAISHIN" in Cartesian coordinates. RAISHIN is based 
on a 3 + 1 formalism of the general relativistic conservation laws of particle number and energy 
momentum, Maxwell's equations, and Ohm's law with no electrical resistance (ideal MHD con- 
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dition) in a curved spacetime (Mizuno et al. 2006). In the RAISHIN code, a conservative, high- 
resolution shock-capturing scheme is employed. The numerical fluxes are calculated using the HLL 
approximate Riemann solver, and flux-interpolated, constrained transport is used to maintain a 
divergence-free magnetic field. The RAISHIN code has proven to be accurate to second order and 
has passed a number of numerical tests including highly relativistic cases and highly magnetized 
cases in both special and general relativity (Mizuno et al. 2006). The RAISHIN code performs 
special relativistic calculations in Minkowski spacetime by changing the metric. 

For our simulations we will choose a force-free helical magnetic field as the initial configuration. 
A force-free configuration is a reasonable choice for the strong magnetic field cases that we study 
here. In general, the force-free equilibrium of a static cylinder is described by the equation 

dB z dB^ 

In particular we choose a poloidal magnetic field component of the form 

R = B ° (?) 

z [1 + {R/a) 2 ] a ' { ' 

for which one finds a toroidal magnetic field component of the form 



B l [l + (R/afY--l-2a(R/af 

* (R/a)[l + (R/a) 2 }«\ 2a - 1 ' { > 

where R is the cylindrical radial position normalized by a simulation scale unit L = 1, Bo parame- 
terizes the magnetic field amplitude, a is the characteristic radius of the column, and a is the pitch 
profile parameter. The pitch profile parameter determines the radial profile of the magnetic pitch 
P = RBz/B^, and provides a measure of the twist of the magnetic field lines. With our choice for 
the force-free field, the magnetic pitch can be written as 



P ~ {R/a? ^ [1 + (R/a) 2 } 2 * - 1 - 2a(R/a) 2 ' (4) 

If the pitch profile parameter a < 1, the magnetic pitch increases with radius. If a > 1, the 
magnetic pitch decreases. When a = 1, the magnetic pitch is constant. This configuration is the 
same as that used in previous non-relativistic work (Appl et al. 2000; Baty 2005). 

The simulation grid is periodic along the axial direction. As a consequence the allowed axial 
wavelengths are restricted to A = L z /n < L z , with n a positive integer and L z is the grid length. 
The grid is a Cartesian (x,y,z) box of size 4L x 4L x L z with grid resolution of AL = L/40. In 
simulations we choose two different column radii, a, relative to the column length, L z : (case A) 
a = (1/16)L 2 = (1/8)L and (case B) a = (1/12)L 2 = (1/4)L. For case AL Z = 2L and for case B 
L z = 3L. In terms of a, the simulation box size is 32a x 32a x 16a for case A and 16a x 16a x 12a 
for case B. We impose outflow boundary conditions on the transverse boundaries at x = y = ±2L 
(x = y = ±16a for case A and x = y = ±8a for case B). 
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We consider a low gas pressure medium with constant p = Po = 0.02poc 2 for the equilibrium 
state, and two different density profiles: (case u) uniform density p = l.Opo and (case n) non- 
uniform density decreasing proportional to the magnetic field strength, p = p\B 2 with p\ = lO.Opo- 
The equation of state is that of an ideal gas with p = (T — l)pe, where e is the specific internal 
energy density and the adiabatic index T = 5/3. The specific enthalpy is h = 1 + e/c 2 + p/pc 2 . The 
magnetic field amplitude is Bq = OA^/ 4ttpqc 2 leading to a low plasma-/? near the axis. The sound 
speed is c s = (Tp/ ph) 1 / 2 and the Alfven speed is va = [B 2 /(ph + B 2 )] 1 / 2 . 

In order to investigate different radial pitch profiles, we perform simulations with: (case 1) 
constant pitch, a = 1, (case 2) increasing pitch, a = 0.75, and (case 3) decreasing pitch a = 2.0. 
The radial profiles of the magnetic field components, the magnetic pitch, and the sound and Alfven 
speeds for the different cases are shown in Figure 1. The radial profile of the magnetic field 
components for cases A and B are the same when normalized by a. 
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Fig. 1. — Radial profile of (a) the toroidal magnetic field (B^),(b) the axial magnetic field (B z ), 
(c) the magnetic pitch, P = RBz/B^, (d) magnetic and gas (dash-dot line) pressure, (e) the rest 
mass density, (f) the sound speed, and (g) the Alfven speed. Uniform density cases are in black 
and declining density cases are in red where (solid) indicates constant pitch, (dotted) indicates 
increasing pitch, and (dashed) indicates decreasing pitch. The radial velocity perturbation profile 
is shown in panel (h). 
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The initial MHD equilibrium configuration is perturbed by a small radial velocity with profile 
given by 

vr = Sv exp (^~~j^J cos(to6>) sin • (5) 

The amplitude of the perturbation is taken to be Sv = 0.01c with radial width R a = 0.5L (4a for 
case A and 2a for case B). We choose m = 1 and n = 1 in the above formula. This is identical to 
imposing (m,n) = (—1,-1), because of the symmetry between (m, n) and (— m, — n) pairs. The 
various different cases that we have considered are listed in Table 1. 
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Table 1. Models and Parameters 



Case 


a/L 


a 


P 


Pitch 


L z /L 


L z /a 


vao/c 


Alu 


0.125 


1.0 


uniform 


constant 


2.0 


16 


0.4 


Aln 


0.125 


1.0 


decrease 


constant 


2.0 


16 


0.3 


Blu 


0.25 


1.0 


uniform 


constant 


3.0 


12 


0.4 


Bin 


0.25 


1.0 


decrease 


constant 


3.0 


12 


0.3 


A2u 


0.125 


0.75 


uniform 


increase 


2.0 


16 


0.4 


A3u 


0.125 


2.0 


uniform 


decrease 


2.0 


16 


0.4 


B2u 


0.25 


0.75 


uniform 


increase 


3.0 


12 


0.4 


B3u 


0.25 


2.0 


uniform 


decrease 


3.0 


12 


0.4 


A2n 


0.125 


0.75 


decrease 


increase 


2.0 


16 


0.3 


A3n 


0.125 


2.0 


decrease 


decrease 


2.0 


16 


0.3 


B2n 


0.25 


0.75 


decrease 


increase 


3.0 


12 


0.3 


B3n 


0.25 


2.0 


decrease 


decrease 


3.0 


12 


0.3 
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Results 



3.1. Time evolution of the CD kink instability 

As an indicator of the growth of the CD kink instability we use the volume-averaged kinetic 
energy transverse to the z-axis within a cylinder of radius R/L < 1.0 written as 



Vh 



P v x + P v y 



dxdydz 



b Jv b 



(6) 



and normalized by the initial volume-averaged transverse magnetic energy. The quantity E\^ n ^ xy 
allows determination of different evolutionary stages (initial exponential (linear growth phase) 
growth, and subsequent non-linear evolution). In all cases, the initial growth regime is characterized 
by an exponential increase in E\.i nxy by several orders of magnitude to a maximum amplitude 
followed by a slow decline in the non-linear regime. The time evolution of the maximum radial 
velocity evaluated over half the grid length, L z /2, shows a growth trend similar to the growth of 
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Fig. 2. — Time evolution of (a) Eki n ^ xy within R/L < 1.0 normalized by the initial volume- 
averaged magnetic energy and (b) the maximum radial velocity normalzed by the Alfven velocity 
on the axis (vao) for constant pitch (a = 1.0). Case A, a = (1/16)L Z : uniform density (solid line) 
and decreasing density (dotted line). Case B, a = (1/12)L Z : uniform density (dashed line) and 
decreasing density (dash-dotted line). Time is in units of = a/ vao, where vao is initial Alfven 
velocity on the axis. Note that normalizing values of E mag)Xy> o, vao and the timescale units, t^, are 
different for uniform and decreasing density cases. 
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Figure 2 shows the time evolution of Ek in ^ xy and maximum radial velocity for cases A and B 
with constant pitch for uniform and decreasing radial density profiles. The effect of changing the 
characteristic radius, a, relative to the grid length, L z , is equivalent to changing the wavelength. 
The wavelength in case A is A = 16 a, and in case B is A = 12 a. Note that according to the Kruskal- 
Shafranov criterion, the instability developes at A > lira. The instability growth rate reaches a 
maximum at X max ~ 10 a, the exact coefficient being dependent on the transverse distribution 
of the density and magnetic pitch. Specifically in the case of constant pitch and uniform density, 
Appl et al. (2000) found A max — 8.43 a and a corresponding growth rate of T max — 0.133 v a o/a. In 
general, one can use the estimate T max 0.1 v a o/a. 

In uniform density (solid and dashed lines in Fig. 2) cases E^i n ^ xy grows exponentially (linear 
growth phase), reaches maximum amplitude at ts ~ lOOt^u for a = (1/16)L 2 and at ts ~ 90£au f° r 
a = (1/12)L 2 , and slowly declines in the non-linear regime. The timescale t is in units of the Alfven 
crossing time of the characteristic radius a, i.e., tAu = a/ vao where the initial Alfven velocity on 
the axis is vao — 0.4 c. The difference in growth is a result of the different wavelengths, longer 
in case A and shorter in case B. Both cases show similar amplitude of Eki n ,xy at transition to the 
non-linear stage. Here we find that the shorter wavelength, A = 12 a, grows slightly more rapidly 
than the longer wavelength, A = 16 a, as would be expected if X m ax < 12 a. 

In decreasing density cases (dotted lines and dash-dotted lines in Fig. 2), the initial growth 
in the linear growth phase is more rapid and reaches maximum amplitude at ts ~ 50tAn f° r 
a = (1/16)L 2 and for a = (1/12)L Z , and slowly declines in the non-linear regime. Here the initial 
Alfven velocity on the axis is vao — 0.3 c and it follows that i^n = (4/3)£au- This more rapid 
growth is a result of dependence of the growth rate on the Alfven velocity. The density decline in 
the decreasing case leads to a more gradual radial decline in the Alfven velocity (see Fig. lg) than 
for the uniform density cases, and on average a higher Alfven velocity. The maximum amplitude 
of Eki n ^ xy relative to the initial volume-averaged magnetic energy for the decreasing density cases 
is almost the same as the uniform density cases. The maximum radial velocities normalized by the 
initial axial Alfven velocity are higher for the decreasing density cases but in absolute terms are 
almost the same. We note that the normalized maximum radial velocity in the decreasing density 
cases appears larger than in the uniform density cases because the normalizing axial Alfven velocity 
is smaller for the decreasing density cases. 

The effect of different radial pitch on the time evolution of E^ in ^ xy and the maximum radial 
velocity for cases with different radial pitch are shown in Figure 3. The three uniform density 
and three decreasing density cases show similar linear and non-linear evolution but with different 
time scales. Growth in the linear stage is more rapid for the decreasing density cases than for 
the uniform density cases. This can be attributed to the higher average Alfven velocity in the 
decreasing density case. For both uniform and decreasing density cases the increasing pitch case 
(a = 0.75) grows more slowly and reaches maximum at a later time with a smaller value for Eki n ,xy 
than the constant pitch case (a = 1.0). On the other hand, the decreasing pitch case (a = 2.0) 
grows more rapidly and reaches maximum at an earlier time with a larger value for Ek in ^ xy than the 
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Fig. 3. — Time evolution of (a,c) E^ ntXy within R/L < 1.0 normalized by the initial volume- 
averaged magnetic energy and (b,d) maximum radial velocity normalzed by the initial Alfven 
velocity on the axis (vao)- Upper panels show the uniform and lower panels show the decreasing 
(lower) density cases with constant pitch a = 1.0 (solid lines), increasing pitch a = 0.75 (dotted 
lines), and decreasing pitch a = 2.0 (dashed lines) cases for a = (1/8)L. Time is in units of a/vAo, 
where vao is initial Alfven velocity on the axis. Note that normalizing values of Emag,xy,o, ^Ao and 
the timescale units, tA, are different for uniform and decreasing density cases. 
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constant pitch case. Although the transition time from linear to non-linear evolution is different 
for each pitch case, the maximum radial velocity is almost the same at transition. The different 
growth rates as a function of the radial pitch profile are consistent with the non-relativistic linear 
analysis in Appl et al. (2000). 
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3.2. Three dimensional structure of the CD kink instability 

The results of the previous subsection show that the instability excites relatively slow motions, 
v <C va, and that the overall energy of the system does not vary too much. Nevertheless, the 
three-dimensional structure of the system is strongly distorted, as will be shown in this subsection. 
We note that results for case B with the shorter wavelength are essentially the same as for case 
A. Moreover, the evolution is similar in the uniform and decreasing density cases and only occurs 
faster in the decreasing density cases. Therefore we show the results of case A with decreasing 
density only. 

Figure 4 shows the time evolution of a density isosurface for the constant pitch case Aln. 
Displacement of the initial force-free helical magnetic field leads to a helically twisted magnetic 
filament around the density isosurface. At longer times, the radial displacement of the high density 
region, best seen in transverse density slices at the grid midplane z = L z /2 (see right panels in 
Fig. 4) , slows significantly. Continuing outwards radial motion is confined to a lower density sheath 
around the high density core. 

The transverse growth of helical twisting is illustrated by the time evolution of magnetic field 
line displacement seen from the pole-on view shown in Figure 5. Here we follow field lines anchored 
initially at R/L = 0.2 at the z = surface. Displacement of the magnetic field rapidly increases 
with time in the linear stage and the displacement continues to gradually increase throughout the 
non-linear phase. Thus, the magnetic field displacement indicates continued growth even though 
the high density region has ceased significant outwards motion by the terminal simulation time. 
Evolution of the helical perturbation in the magnetic field is accompanied by an appropriate evolu- 
tion of the current distribution. Initially a large axial current is located near the axis. In the linear 
growth phase, the axial current is displaced and twists helically. The axial current twist is cospatial 
with the high density twist. In the linear growth phase the axial current gradually decreases and 
in the non-linear phase remains constant in magnitude. 

Figure 6 shows density isosurfaces and transverse density slices at the grid midplane, z = L z /2, 
for the cases A2n (increasing pitch) and A3n (decreasing pitch) that can be compared to Figure 4 
for case Aln with constant pitch. The 3D density structure from A2n looks similar to that from 
Aln. The CD kink instability grows exponentially initially. The transverse density slice at the 
grid midplane (Fig. 6c) shows little outwards motion of the high density region in the non-linear 
stage similar to the constant pitch case Aln, however, there is little outwards motion in the low 
density sheath surrounding the high density region. This result is somewhat different from the 
constant pitch case Aln and suggests a significant reduction in kink amplitude growth. Results 
from the decreasing pitch case A3n are very different. Figure 6b shows a more slender helical 
density structure wrapped by the magnetic field. Recall that transition from linear to non-linear 
evolution was reached at a much earlier time, ts ~ 50t^n ( see Fig. 3), than is shown here. While 
the density cross section (Fig. 6d) is similar to that of the constant pitch case Aln (Fig. 4f), radial 
motion continues in the non-linear stage to the end of the simulation. 




Fig. 4. — 3D density isosurface with transverse slices at z = (left) and 2D transverse slices at 
the grid midplane z = L z /2 (right) for case Aln. Colors show the logarithm of the density with 
magnetic field (white) lines. 



-14- 



Case: Aln (Non-uniform Denisity, Constant Pitch, a=l/8L) 
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Fig. 5. — Magnetic field structure seen from a pole on view for case Aln. 




Fig. 6. — 3D density isosurface with transverse slices at z = (upper panels) and 2D transverse 
slices at the grid midplane z = L z /2 (lower panels) for case A2n (left) and A3n (right). Colors 
show the logarithm of the density with magnetic field (white) lines. 




Fig. 7. — 3D axial current, j z , isosurface (upper panels) with magnetic field (white) lines and 
(bottom panels) magnetic field structure seen from a pole on view for cases A2n (left) and A3n 
(right). 
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Figure 7 shows the axial current isosurface and the magnetic field line displacement seen from 
the pole-on view for cases A2n and A3n. Like the constant pitch decreasing density case Aln, the 
axial current twist is cospatial with the high density twist. In the non-linear stage, the helically 
twisted axial current does not change significantly in case A2n. A pole on view of field lines 
anchored initially at R/L = 0.2 at the z = surface (see Fig. 7c) shows that the magnetic field 
displacement does not increase significantly in the non-linear stage. This along with the density 
results shown in Figure 6a & c suggest significant non-linear stabilization of the kink mode in the 
presence of increasing pitch. On the other hand, the axial current isosurface for the decreasing 
pitch case A3n seen in Figure 7b makes a slender helical structure coincident with the high density 
structure. Large displacement in the magnetic field can be seen in the pole on view (Fig. 7d). 



4. Summary and Discussion 

We have investigated development of the CD kink instability of force-free helical magnetic 
equilibria in order to study the influence of the initial configuration on the linear and nonlinear 
evolution. We followed the time evolution where the plasma column was perturbed to excite the 
m = ±1 azimuthal mode. The kink developed initially as predicted by linear stability analysis. 
The rate at which the kink developed and, in particular, the non-linear behavior depended on the 
density and magnetic pitch radial profile. For a constant magnetic pitch profile a density decline 
with radius led to faster linear growth and made a transition to the non-linear stage sooner than 
occured for a uniform density as a result of the more gradual decrease in the Alfven velocity with 
radius. While there were some differences in the kink density structure in these two different cases, 
amplitude growth of the kink continued to the terminal simulation time in both cases. 

More profound differences in development accompanied different magnetic pitch profiles. An 
increase in the magnetic pitch with radius led to a reduced growth rate in the linear (exponential) 
growth stage and made a transition to non-linear behavior at a later time when compared to the 
cases with constant magnetic pitch. In the nonlinear stage, amplitude growth of the kink appeared 
to have ceased by the terminal simulation time. On the other hand, a decrease in the magnetic 
pitch led to more rapid growth in the the linear growth stage and made a transition to non-linear 
behavior at an earlier time when compared to the cases with constant magnetic pitch. In the non- 
linear stage, a slender helical density and current carrying plasma column wrapped by magnetic 
field developed and amplitude growth continued to increase throughout the simulation. Thus, the 
linear growth and nonlinear evolution of the CD kink instability depended on the radial density 
profile and strongly depended on the magnetic pitch profile. Increasing magnetic pitch with radius 
was clearly stabilizing. 

In all cases, the initial axisymmetric structure is strongly distorted by the kink instability, 
even though not disrupted. It is important to note that the instability develops relatively slowly, 
that is, the characteristic time for the instability to affect strongly the initial structure is roughly 
r ~ 100 va/cl- The growth rate of the kink instability is roughly k ~ 0.1 va/c, the exact coefficient 
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being dependent on the structure of the undisturbed state. Therefore it is not accidental that the 
full development of the instability takes a long time. In a jet context, our simulations correspond 
to a perturbation that remains at rest in the flow frame. Therefore, in order to check whether 
the instability would affect a jet flow, one has to compare r with the propagation time. In the 
relativistically moving case, time dilation slows the instability further by the jet Lorentz factor so 
that the condition for the instability to affect the jet structure may be written as 



In relativistic jets, one could take va ~ c. In order to find the final criterion, one has to know 
how the jet radius, a, and Lorentz factor, 7, vary with the distance, z. If the jet is narrow enough 
so that Qa 2 /c < z, where Q is the angular velocity at the base of the jet, one can use the scaling 
(Tschekhovskoy et al. 2008; Komissarov et al. 2009; Lyubarsky 2009) 7 ~ Qa. In this case one 
finds that the criterion for the kink instability is 



This means that the instability could affect the jet structure only if the jet expands slowly enough. 
Assuming the parabolic shape for the jet, Qa/c = £(Qz/c) k , where k < 1 and £ ~ 1 are dimensionless 
numbers, one finds that the instability develops only if k < 1/2. Even in this case, the characteristic 
scale for the development of the instability is very large, 



The above estimates should be considered as preliminary because we assumed that the jet 
moves as a whole so that there is a frame of reference in which the plasma is at rest. Clearly, 
the criterion for the CD kink instability can be modified by the effects of relativistic rotation, 
gradual shear, a surrounding sheath, and sideways expansion. For example, in relativistic Poynting 
dominated jets the growth rate of the instability strongly depends on the transverse profile of the 
poloidal magnetic field so that r goes to infinity when the poloidal field is uniform (Istomin & 
Pariev 1994, 1996; Lyubarskii 1999). On the other hand, Poynting dominated jets are accelerated 
such that the acceleration zone spans a large range of scales (Tschekhovskoy et al. 2008; Komissarov 
et al. 2008; Lyubarsky 2009). Therefore, a small kink growth rate does not necessarily preclude 
significant growth of helical perturbations. We plan to perform 3D RMHD simulations including 
some of the missing physical effects such as jet flow and rotation in future work. In particular, this 
will allow us to investigate the interaction between CD and KH driven structure in the relativistic 
regime. 
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NASA NNG05GK73G, NNX07AJ88G, and NNX08AG83G. This work is supported in part by 
US-Israeli Binational Science Foundation grant 2006170. The simulations were performed on the 
Columbia Supercomputer at NAS Division at NASA Ames Research Center and the Altix3700 BX2 
at YITP in Kyoto University. 
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Numerical Tests 



In the simulations presented here numerical viscosity and resistivity are associated with the 
finite grid resolution. It is generally expected that numerical viscosity affects the quantitative 
growth or damping of instabilities and numerical resistivity allows unexpected magnetic reconnec- 
tion even when we assume the ideal MHD condition. Therefore we perform several numerical tests 
to check whether saturation is adequately captured and investigate how the results in the non-linear 
evolution stage depend on numerical resolution. 

Case: Alu (Uniform density Constant pitch, a=l/8L) 
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Fig. 8. — Time evolution of (a) Eki n ,xy normalized by initial volume-averaged magnetic energy 
within R/L < 1.0 and (b) maximum radial velocity normalzed by initial Alfven velocity on the axis 
for case Alu with grid resolution of AL = L/60 (dashed), L/40 (solid), L/30 (dotted), and L/20 
(dash-dotted). 



First, we repeat case Alu (uniform density with constant pitch) for four grid resolutions from 
20 to 60 computational zones per simulation length unit L = 8a. Figure 8 shows the time evolution 
of Ekin )Xy and the maximum radial velocity for different grid resolutions. The results depend on 
grid resolution. Higher grid resolution shows faster growth and an earlier transition time to the 
non-linear stage at a larger amplitude. This result is the same as that found in the resolution 
study of Baty & Keppens (2002). The lowest grid resolution with AL = L/20 is clearly inadequate. 
Although growth does depend on grid resolution, the difference between our choice of AL = L/40 
and the highest grid resolution of AL = L/60 is not significant. Therefore our simulation results 
adequately capture the CD kink instability. 

Second, we check the influence of the transverse boundary. Figure 9 shows the time evolution 
of Ekin,xy and the maximum radial velocity within a box of length L 2 /2 for a simulation box with 
transverse boundaries at x = y = ±2L (±16a) and a simulation box with transverse boundaries 



Case: Alu (Uniform density Constant pitch, a=l/8L) 
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Fig. 9. — Time evolution of (a) E^^y normalized by initial volume-averaged magnetic energy 
within R/L < 1.0 and (b) maximum radial velocity normalzed by initial Alfven velocity on the axis 
for case Alu with transverse boundaries at x = y = ±2L (±16a: solid) and x = y = ±3L (±24a: 
dashed). 



at x = y = ±3L (±24a). The results for the time evolution of E^^y and the maximum radial 
velocity are the same in both cases. We conclude that behavior of the CD kink instability is not 
affected by our choice for the transverse boundary. 

Third, we check the influence of the initial gas pressure. In the force-free model we ignore the 
inertia and pressure of the plasma. However, in the simulation the relativistic ideal MHD equations 
are solved numerically not the force-free equations. Figure 10 shows the time evolution of E^ n ^ xy 
and the maximum radial velocity in a box of length L z /2 for an initial pressure of po = 0.02poc 2 
and po = 0.005/)oc 2 . The results of the time evolution of E^i n ^ xy and the maximum radial velocity 
show the same linear growth rate in both cases. However, the lower initial pressure case reaches 
saturation at a slightly later time with higher amplitude but with a larger decline in the non-linear 
stage. Thus, our required non-zero initial gas pressure does affect the saturation level and non-linear 
development of the CD kink instability somewhat. 
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